/*
* INSANE - Interactive Structural Analysis Environment
*
* Copyright (C) 2003-2005
* Universidade Federal de Minas Gerais
* Escola de Engenharia
* Departamento de Engenharia de Estruturas
* 
* Author's email :    insane@dees.ufmg.br
* Author's website :  http://www.dees.ufmg.br/insane
* 
* This program is free software; you can redistribute it and/or
* modify it under the terms of the GNU General Public License
* as published by the Free Software Foundation; either version 2
* of the License, or any later version.
* 
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
* GNU General Public License for more details.
* 
* You should have received a copy of the GNU General Public License
* along with this program; if not, write to the Free Software
* Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
*/
package br.ufmg.dees.insane.model.femModel.problemDriver;

import java.io.Serializable;
import java.util.ArrayList;
import java.util.ListIterator;

import br.ufmg.dees.insane.analysisModel.Beam;
import br.ufmg.dees.insane.analysisModel.FrameElmAnalysis;
import br.ufmg.dees.insane.analysisModel.Grid;
import br.ufmg.dees.insane.analysisModel.PlaneFrame;
import br.ufmg.dees.insane.analysisModel.PlaneTruss;
import br.ufmg.dees.insane.analysisModel.SpaceFrame;
import br.ufmg.dees.insane.analysisModel.SpaceTruss;
import br.ufmg.dees.insane.model.femModel.element.Element;
import br.ufmg.dees.insane.model.femModel.element.FrameElement;
import br.ufmg.dees.insane.model.femModel.element.force.ElementForce;
import br.ufmg.dees.insane.model.femModel.node.Node;
import br.ufmg.dees.insane.model.femModel.pointforce.PointForce;
import br.ufmg.dees.insane.shape.Shape1DCart;
import br.ufmg.dees.insane.util.IMatrix;
import br.ufmg.dees.insane.util.IPoint3d;
import br.ufmg.dees.insane.util.IVector;

/**
 * Class representing FrameProblemDriver
 * @author Samir, Flavio, Roque
 * @since 20/09/2006
 */
public class FrameProblemDriver implements ProblemDriver, Serializable {

	private static final long serialVersionUID = 1L;

	/* (non-Javadoc)
     * @see br.ufmg.dees.insane.model.femModel.problemDriver.ProblemDriver#getC(br.ufmg.dees.insane.model.femModel.element.Element)
     */
    public IMatrix getC(Element e) {
    	
    	FrameElement e1 = (FrameElement)e;
    	
		//Cartesian coordinates matrix
    	IPoint3d point1 = ((Node)e1.getIncidence().get(0)).getPoint();
		IPoint3d point2 = ((Node)e1.getIncidence().get(1)).getPoint();
		IMatrix cn = new IMatrix(3,3);
    	cn.setElement(0, 0, point1.getX());
    	cn.setElement(0, 1, point1.getY());
    	cn.setElement(0, 2, point1.getZ());
    	cn.setElement(1, 0, point2.getX());
    	cn.setElement(1, 1, point2.getY());
    	cn.setElement(1, 2, point2.getZ());
    	
    	//Nodal angles
    	double ang1 = ((Node)e1.getIncidence().get(0)).getAngle();
    	double ang2 = ((Node)e1.getIncidence().get(1)).getAngle();
    	
    	IMatrix t = e1.getAnalysisModel().getTransformationMatrix(null, cn);
    	IMatrix k = this.localCorrectedStiffnessMatrix(e1);
		IMatrix r = ((FrameElmAnalysis)e1.getAnalysisModel()).nodalRotationMatrix(ang1, ang2);
		IMatrix aux1 = new IMatrix(r.getNumCol(), t.getNumCol());
		aux1.mulTransposeRight(r, t);
		IMatrix aux2 = new IMatrix(aux1.getNumCol(), k.getNumRow());
		aux2.mul(aux1, k);
		aux1.mul(aux2, t);
		aux2.mulTransposeRight(aux1, r);
		
		return (aux2);
		
    }

    /* (non-Javadoc)
     * @see br.ufmg.dees.insane.model.femModel.problemDriver.ProblemDriver#getR(br.ufmg.dees.insane.model.femModel.element.Element)
     */
	public IVector getE(Element e) {
	    IVector E = new IVector(((FrameElement)e).getNumberOfDegreesOfFreedom());
	    return E;
	    
	}

    /* (non-Javadoc)
     * @see br.ufmg.dees.insane.model.femModel.problemDriver.ProblemDriver#getF(br.ufmg.dees.insane.model.femModel.element.Element)
     */
    public IVector getF(Element e){
    	
    	FrameElement e1 = (FrameElement)e;
    	FrameElmAnalysis anl = (FrameElmAnalysis)e1.getAnalysisModel();
    	
    	//Cartesian coordinates matrix
    	IPoint3d point1 = ((Node)e1.getIncidence().get(0)).getPoint();
		IPoint3d point2 = ((Node)e1.getIncidence().get(1)).getPoint();
		IMatrix cn = new IMatrix(3,3);
    	cn.setElement(0, 0, point1.getX());
    	cn.setElement(0, 1, point1.getY());
    	cn.setElement(0, 2, point1.getZ());
    	cn.setElement(1, 0, point2.getX());
    	cn.setElement(1, 1, point2.getY());
    	cn.setElement(1, 2, point2.getZ());
    	
    	//Nodal angles
    	double ang1 = ((Node)e1.getIncidence().get(0)).getAngle();
    	double ang2 = ((Node)e1.getIncidence().get(1)).getAngle();
    	
		IMatrix t = anl.getTransformationMatrix(null, cn);
		IVector f = this.localCorrectedEquivalentForceVector(e1);
		IMatrix r = anl.nodalRotationMatrix(ang1, ang2);
		IVector aux = new IVector(f.getSize());
		t.transpose();
		aux.mul(t, f);
		IVector aux1 = new IVector(f.getSize());
		aux1.mul(r, aux);
		return(aux1);
		
    }

    /* (non-Javadoc)
     * @see br.ufmg.dees.insane.model.femModel.problemDriver.ProblemDriver#getStateVariable()
     */
    public String getStateVariableLabel() {
      return "Displacements";
    }

    /* (non-Javadoc)
     * @see br.ufmg.dees.insane.model.femModel.problemDriver.ProblemDriver#getDualVariable()
     */
    public String getDualVariableLabel() {
      return "Loads";
    }

    /* (non-Javadoc)
     * @see br.ufmg.dees.insane.model.femModel.problemDriver.ProblemDriver#getInternalVariable()
     */
    public String getInternalVariableLabel() {
      return "Epsilon";
    }

    /* (non-Javadoc)
     * @see br.ufmg.dees.insane.model.femModel.problemDriver.ProblemDriver#getInternalDualVariable()
     */
    public String getInternalDualVariableLabel() {
      return "Sigma";
    }
	
    /* (non-Javadoc)
     * @see br.ufmg.dees.insane.model.femModel.problemDriver.ProblemDriver#getMPInternalVariableLabel()
     */
	public String getMPInternalVariableLabel() {
		return "Strains";
	}

    /* (non-Javadoc)
     * @see br.ufmg.dees.insane.model.femModel.problemDriver.ProblemDriver#getMPInternalDualVariableLabel()
     */
	public String getMPInternalDualVariableLabel() {
		return "Stress";
	}
    
	/** Returns the corrected local stiffness matrix of the specified FrameElement, accordingly to this AnalysisModel type.
	*This matrix is corrected to consider the liberations at the extremities of the FrameElement.
	*@param elm The FrameElement whose corrected local stiffness matrix is desired.
	*@return The corrected local stiffness matrix of the specified FrameElement, accordingly to this AnalysisModel type.
	*/
	public IMatrix localCorrectedStiffnessMatrix(FrameElement elm) {
		int ndf = elm.getNumberOfDegreesOfFreedom();
		FrameElmAnalysis anl = (FrameElmAnalysis)elm.getAnalysisModel();
		
		IMatrix mat = new IMatrix(ndf, ndf);
		IMatrix mat2 = new IMatrix(ndf, ndf);
		mat2.setZero();
		
		double[] prop = elm.getPropertiesArray();
		
		mat = anl.getLocalStiffnessMatrix(prop);
		
		//Changing stiffness matrix because of liberations at extremity
		int l=-1;
		double kmj, klj, kll, kml;
		boolean[] liberations = elm.getLiberations();
		for (int i=0; i<ndf; i++) {
			l++;
			if (liberations[i]) {
				kll = mat.getElement(l, l);
				for (int j=0; j<ndf; j++) {
					for (int m=0; m<ndf; m++) {
						kmj = mat.getElement(m, j);
						klj = mat.getElement(l ,j);
						kml = mat.getElement(m, l);
						mat2.setElement(m, j, kmj - (klj / kll) * kml);
					}
				}
				mat.set(mat2);
			}
		}
		
		return (mat);
		
	}
	
	/** Returns the local equivalent nodal force vector of the specified FrameElement, accordingly to this AnalysisModel type.
	*This vector does not consider the possible liberations at the extremities of the FrameElement.
	*@param elm The FrameElement whose local equivalent nodal force vector is desired.
	*@return The local equivalent nodal force vector of the specified FrameElement, accordingly to this AnalysisModel type.
	*/
	public IVector localEquivalentForceVector(FrameElement elm) {
		int ndf = elm.getNumberOfDegreesOfFreedom();
		int ndfByNode = elm.getAnalysisModel().getDOFLabels(1).size();
		FrameElmAnalysis anl = (FrameElmAnalysis)elm.getAnalysisModel();
		Shape1DCart shp = (Shape1DCart)elm.getShape();
		
		//Cartesian coordinates matrix
		IPoint3d point1 = ((Node)elm.getIncidence().get(0)).getPoint();
		IPoint3d point2 = ((Node)elm.getIncidence().get(1)).getPoint();
		IMatrix cn = new IMatrix(3,3);
    	cn.setElement(0, 0, point1.getX());
    	cn.setElement(0, 1, point1.getY());
    	cn.setElement(0, 2, point1.getZ());
    	cn.setElement(1, 0, point2.getX());
    	cn.setElement(1, 1, point2.getY());
    	cn.setElement(1, 2, point2.getZ());
    	
		IVector lEqForce = new IVector(ndf);
		IVector a = new IVector(ndfByNode);
		IVector b = new IVector(ndfByNode);
		IVector v1;
		IVector v2;
		IMatrix v3;
		IVector f1 = new IVector(ndf);
		IVector f2 = new IVector(ndf);
		IMatrix n;
		IMatrix n1;
		IMatrix n2;
		
		//Calculating equivalent forces for elmPointForces
		PointForce elmPF;
		ListIterator elmPFs = ((ArrayList)elm.getElementValues().get(FrameElement.ELM_POINT_FORCES)).listIterator();
		while (elmPFs.hasNext()) {
			elmPF = (PointForce) elmPFs.next();
			
			//Setting the values for vector a: (fe = N * a)
			for (int i=0; i<ndfByNode; i++) {
				double scale = elmPF.getLoadCase().getFactor();
				a.setElement(i, elmPF.getForce(i)*scale);
			}
			
			v1 = shp.getShapeFunction(null, cn);
			v3 = shp.getDerivedShapeFunction(null, cn);
			n = anl.getStateVariableOperator(v3, v1, null);
			n.transpose();
			f1.mul(n, a);
			lEqForce.add(f1);
		}
		
		//Calculating equivalent forces for elementForces
		a.zero();
		ElementForce bf;
		PointForce pf1;
		PointForce pf2;
		ListIterator bfs = ((ArrayList)elm.getElementValues().get(FrameElement.ELM_LINE_FORCES)).listIterator();
		while (bfs.hasNext()) {
			bf = (ElementForce) bfs.next();
			pf1 = (PointForce) bf.getPointForcesList().get(0);
			pf2 = (PointForce) bf.getPointForcesList().get(1);
			
			//Setting the values of the vectors a and b (lEFV = IntegratedShapeFunctionMultX * a + IntegratedShapeFunction * b)
			for (int i=0; i<ndfByNode; i++) {
				double scale = pf1.getLoadCase().getFactor();
				a.setElement(i, (pf2.getForce(i) - pf1.getForce(i))*scale / (pf2.getCoord().x - pf1.getCoord().x));
				b.setElement(i, -(a.getElement(i)) * pf2.getCoord().x + pf2.getForce(i));
			}
			
			//Setting the arrays which contain the coordinates of the integration limits
			double[] a1 = new double[3];
			a1[0] = pf1.getCoord().getX();
			a1[1] = pf1.getCoord().getY();
			a1[2] = pf1.getCoord().getZ();
			double[] a2 = new double[3];
			a2[0] = pf2.getCoord().getX();
			a2[1] = pf2.getCoord().getY();
			a2[2] = pf2.getCoord().getZ();
			
			v1 = shp.getIntegratedShapeFunctionTimesX(a1, a2, cn);
			v2 = shp.getIntegratedShapeFunction(a1, a2, cn);
			
			n1 = anl.getStateVariableOperator(null, v1, null);
			n2 = anl.getStateVariableOperator(null, v2, null);
			n1.transpose();
			n2.transpose();
			f1.mul(n1, a);
			f2.mul(n2, b);
			f1.add(f2);
			lEqForce.add(f1);
		}
		
		//Calculating equivalent forces for temperature changes
		if (elm.getHeight()!=0) {
			f1.zero();
			f1 = new IVector(ndf);
			//double[] deltaTemp = elm.getElementValues(FrameElement.DELTA_TEMP).getDoublePointValue(0);
			double factor1 = elm.getElasticity() * elm.getThermalCoefficient() * elm.getArea() * elm.getElementValues(FrameElement.DELTA_TEMP).getDoublePointValue(1);
			double factor2 = elm.getElasticity() * elm.getThermalCoefficient() * elm.getIz() / elm.getHeight() * (elm.getElementValues(FrameElement.DELTA_TEMP).getDoublePointValue(2) - elm.getElementValues(FrameElement.DELTA_TEMP).getDoublePointValue(0));
			if (anl instanceof SpaceFrame) {
				f1.setElement(0, -factor1);
				f1.setElement(5, -factor2);
				f1.setElement(6, factor1);
				f1.setElement(11, factor2);
			}
			else if (anl instanceof SpaceTruss) {
				f1.setElement(0, -factor1);
				f1.setElement(3, factor2);
			}
			else if (anl instanceof Grid) {
				f1.setElement(2, -factor2);
				f1.setElement(5, factor2);
			}
			else if (anl instanceof Beam) {
				f1.setElement(1, -factor2);
				f1.setElement(3, factor2);
			}
			else if (anl instanceof PlaneFrame) {
				f1.setElement(0, -factor1);
				f1.setElement(2, -factor2);
				f1.setElement(3, factor1);
				f1.setElement(5, factor2);
			}
			else if (anl instanceof PlaneTruss) {
				f1.setElement(0, -factor1);
				f1.setElement(3, factor1);
			}
			lEqForce.add(f1);
		}
		
		//Calculating equivalent forces for prescribed displacements on elements
		f1.zero();
		n1 = localCorrectedStiffnessMatrix(elm);
		v1 = elm.getElmPreDisplacement();
		f2.zero();
		for (int i=0; i<ndfByNode; i++) {
				f2.setElement(i, v1.getElement(i));
				f2.setElement(i+(ndfByNode/2), v1.getElement(i+ndfByNode));
		}
		f1.mul(n1, f2);
		lEqForce.add(f1);
		
		//Adding initial nodal force
		f1 = new IVector(ndf);
		double[] iniEqforces = elm.getInitialEquivalentNodalForce();
		for (int i=0; i<ndfByNode; i++) {
			f1.setElement(i, iniEqforces[i]);
			f1.setElement(i+ndfByNode, iniEqforces[i+ndfByNode]);
		}
		
		lEqForce.add(f1);
		
		return (lEqForce);
		
	}
	
	/** Returns the corrected local equivalent nodal force vector of the specified FrameElement, accordingly to this AnalysisModel type.
	*This matrix is corrected to consider the liberations at the extremities of the FrameElement.
	*@param elm The FrameElement whose corrected local equivalent nodal force vector is desired.
	*@return The corrected local equivalent nodal force vector of the specified FrameElement, accordingly to this AnalysisModel type.
	*/
	public IVector localCorrectedEquivalentForceVector(FrameElement elm) {
		int ndf = elm.getAnalysisModel().getDOFLabels(2).size();
		FrameElmAnalysis anl = (FrameElmAnalysis)elm.getAnalysisModel();
		
		IVector vec = new IVector(ndf);
		vec = this.localEquivalentForceVector(elm);
		
		IVector vec2 = new IVector(ndf);
		
		//FrameElement properties
		double[] prop = elm.getPropertiesArray();
		
		IMatrix mat = new IMatrix (ndf, ndf);
		mat = anl.getLocalStiffnessMatrix(prop);
		
		IMatrix mat2 = new IMatrix(ndf, ndf);
		mat2.setZero();
		
		int l=-1;
		double kmj, klj, kll, kml, fl, fm;
		boolean[] liberations = elm.getLiberations();
		for (int i=0; i<ndf; i++) {
			l++;
			if (liberations[i]) {
				kll = mat.getElement(l, l);
				fl = vec.getElement(l);
				
				for (int j=0; j<ndf; j++) {
					for (int m=0; m<ndf; m++) {
						kmj = mat.getElement(m, j);
						klj = mat.getElement(l ,j);
						kml = mat.getElement(m, l);
						fm = vec.getElement(m);
						vec2.setElement(m, kml / kll * fl - fm);
						mat2.setElement(m, j, kmj - (klj / kll) * kml);
					}
				}
				vec2.negate();
				vec.set(vec2);
				mat.set(mat2);
			}
		}
		return (vec);	
	}
	
	public IMatrix getA(Element e) {
		// TODO Auto-generated method stub
		return null;
	}

	public IMatrix getB(Element e) {
		// TODO Auto-generated method stub
		return null;
	}

	public IMatrix getIncrementalA(Element e) {
		// TODO Auto-generated method stub
		return null;
	}

	public IMatrix getIncrementalB(Element e) {
		// TODO Auto-generated method stub
		return null;
	}

	public IMatrix getIncrementalC(Element e) {
		// TODO Auto-generated method stub
		return null;
	}

	public IMatrix getTotalA(Element e) {
		// TODO Auto-generated method stub
		return null;
	}

	public IMatrix getTotalB(Element e) {
		// TODO Auto-generated method stub
		return null;
	}

	public IMatrix getTotalC(Element e) {
		// TODO Auto-generated method stub
		return null;
	}

	public String getUnknownDualVariableLabel() {
		// TODO Auto-generated method stub
		return null;
	}
	
}
